#include "dyn_ode.h"

ODE::ODE(MultiFunctor &func,
         colvec _initCond,
         double _start,
         double _end,
         double _initStepsize,
         double _minStepsize,
         double _maxStepsize,
         double _absTol,
         double _relTol,
         int _maxIterations)
    : integrand(func)
{
    setInitCond(_initCond);
    setStart(_start);
    setEnd(_end);
    setInitStepsize(_initStepsize);
    setMinStepsize(_minStepsize);
    setMaxStepsize(_maxStepsize);
    setAbsTol(_absTol);
    setRelTol(_relTol);
    setMaxIterations(_maxIterations);
}

int ODE::integrate()
{
    //cout << "Preparing integration...";
    // We determine the dimension of the system
    // as well as the number of steps we are going to take.
    int N = (int) fabs((end-start)/initStepsize)+1;
    int dimension = initCond.n_rows;

    // We prepare a matrix of the proper size.
    // The first column contains the independent variable.
    mat tempResults = mat(N+1,dimension+1);

    // The first row contains the initial data.
    tempResults(0,0) = start;
    tempResults.submat(span(0,0),span(1,dimension)) = initCond.t();

    // Now we compute the complete integral.
    //cout << "Integrating." << endl;
    for (int i=1;i<=N;i++)
    {
        tempResults(i,0) = tempResults(i-1,0)+initStepsize;
        tempResults.submat(span(i,i),span(1,dimension)) = next(tempResults.submat(span(i-1,i-1),span(1,dimension)).t()).t();
    }

    // We store the results in the
    // member variable.
    results = tempResults;

    return 0;
}
